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Abstract 

We consider stochastically modeled chemical reaction systems with mass-action kinetics 
and prove that a product-form stationary distribution exists for each closed, irreducible subset 
of the state space if an analogous deterministically modeled system with mass-action kinet- 
ics admits a complex balanced equilibrium. Feinberg's deficiency zero theorem then implies 
that such a distribution exists so long as the corresponding chemical network is weakly re- 
versible and has a deficiency of zero. The main parameter of the stationary distribution for the 
stochastically modeled system is a complex balanced equilibrium value for the corresponding 
deterministically modeled system. We also generalize our main result to some non-mass-action 
kinetics. 



1 Introduction 

There are two commonly used models for chemical reaction systems: discrete stochastic models 
in which the state of the system is a vector giving the number of each molecular species, and 
continuous deterministic models in which the state of the system is a vector giving the concentra- 
tion of each molecular species. Discrete stochastic models are typically used when the number of 
molecules of each chemical species is low and the randomness inherent in the making and breaking 
of chemical bonds is important. Conversely, deterministic models are used when there are large 
numbers of molecules for each species and the behavior of the concentration of each species is 
well approximated by a coupled set of ordinary differential equations. 

Typically, the goal in the study of discrete stochastic systems is to either understand the evolu- 
tion of the distribution of the state of the system or to find the long term stationary distribution of the 
system, which is the stochastic analog of an equilibrium point. The Kolmogorov forward equation 
(chemical master equation in the chemistry literature) describes the evolution of the distribution 
and so work has been done in trying to analyze or solve the forward equation for certain classes 
of systems ([20]). However, it is typically an extremely difficult task to solve or even numeri- 
cally compute the solution to the forward equation for all but the simplest of systems. Therefore, 
simulation methods have been developed that will generate sample paths so as to approximate the 
distribution of the state via Monte Carlo methods. These simulation methods include algorithms 
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that generate statistically exact (E|22l|23l|2l]|) and approximate ([H SIM 13) sample paths. On 
the other hand, the continuous deterministic models, and in particular mass-action systems with 
complex balancing states, have been analyzed extensively in the mathematical chemistry litera- 
ture, starting with the works of Horn, Jackson, and Feinberg ( 11261 l27l l28l [EH), and continuing 
with Feinberg's deficiency theory in ([161 111 EEU [HI)- Such models have a wide range of appli- 
cations in the physical sciences, and now they are beginning to play an important role in systems 
biology ( |fT3l [231 l37ll ). Recent mathematical analysis of continuous deterministic models has fo- 
cused on their potential to admit multiple equilibria ( liTTl[T2"1 ) and on dynamical properties such as 
persistence and global stability (llJ7ll7ll2~ll4ll6ln. 

One of the major theorems pertaining to deterministic models of chemical systems is the defi- 
ciency zero theorem of Feinberg ( IfTTllTo'lO . The deficiency zero theorem states that if the network 
of a system satisfies certain easily checked properties, then within each compatibility class (invari- 
ant manifold in which a solution is confined) there is precisely one equilibrium with strictly positive 
components, and that equilibrium is locally asymptotically stable ( ifTTlfToTO . The surprising aspect 
of the deficiency zero theorem is that the assumptions of the theorem are completely related to the 
network of the system whereas the conclusions of the theorem are related to the dynamical proper- 
ties of the system. We will show in this paper that if the conditions of the deficiency zero theorem 
hold on the network of a stochastically modeled chemical system with quite general kinetics, then 
there exists a product-form stationary distribution for each closed, irreducible subset of the state 
space. In fact, we will show a stronger result: that a product-form stationary distribution exists so 
long as there exists a complex balanced equilibrium for the associated deterministically modeled 
system. However, the equilibrium values guaranteed to exist by the deficiency zero theorem are 
complex balanced and so the conditions of that theorem are sufficient to guarantee the existence 
of the product-form distribution. Finally, the main parameter of the stationary distribution will be 
shown to be a complex balanced equilibrium value of the deterministically modeled system. 

Product-form stationary distributions play a central role in the theory of queueing networks 
where the product-form property holds for a large, naturally occurring class of models called Jack- 
son networks (see, for example, 11301 , Chapter 3, and IfTOl , Chapter 2) and a much larger class of 
quasi-reversible networks ((201, Chapter 3, OH, Chapter 4, [36|, Chapter 8). Kelly, [30], Section 
8.5, recognizes the possible existence of product-form stationary distributions for a subclass of 
chemical reaction models and gives a condition for that existence. That condition is essentially the 
complex balance condition described below, and our main result asserts that for any mass-action 
chemical reaction model the conditions of the deficiency zero theorem ensure that this condition 
holds. 

The outline of the paper is as follows. In Section [2] we formally introduce chemical reaction 
networks. In Section [3] we develop both the stochastic and deterministic models of chemical re- 
action systems. Also in Section [3] we state the deficiency zero theorem for deterministic systems 
and present two theorems that are used in its proof and that will be of use to us. In Section |4] 
we present the first of our main results: that every closed, irreducible subset of the state space of a 
stochastically modeled system with mass-action kinetics has a product-form stationary distribution 
if the chemical network is weakly reversible and has a deficiency of zero. In Section[5]we present 
some examples of the use of this result. In Section [6] we extend our main result to systems with 
more general kinetics. 
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2 Chemical reaction networks 



Consider a system with m chemical species, {Si, . . . , S m }, undergoing a finite series of chemical 
reactions. For the A;th reaction, denote by v k , i/ k E Z> the vectors representing the number of 
molecules of each species consumed and created in one instance of that reaction, respectively. We 
note that if v k — then the A;th reaction represents an input to the system, and if v' k — then it 
represents an output. Using a slight abuse of notation, we associate each such v k (and u' k ) with 
a linear combination of the species in which the coefficient of Si is v ik , the ith element of v k . 
For example, if v k = [1, 2, 3] T for a system consisting of three species, we associate with v k 
the linear combination S\ + 2S 2 + 3S 3 . For v k = 0, we simply associate u k with 0. Under this 
association, each v k (and u' k ) is termed a complex of the system. We denote any reaction by the 
notation v k — > v' k , where v k is the source, or reactant, complex and v' k is the product complex. We 
note that each complex may appear as both a source complex and a product complex in the system. 
The set of all complexes will be denoted by {u k } := Li k ({u k } U {v' k }). 

Definition 2.1. Let S = {Si}, C = {v k }, and 1Z = {v k — > u' k } denote the sets of species, 
complexes, and reactions, respectively. The triple {S, C, 1Z] is called a chemical reaction network. 

The structure of chemical reaction networks plays a central role in both the study of stochasti- 
cally and deterministically modeled systems. As alluded to in the Introduction, it will be conditions 
on the network of a system that guarantee certain dynamical properties for both models. Therefore, 
the remainder of this section consists of definitions related to chemical networks that will be used 
throughout the paper. 

Definition 2.2. A chemical reaction network, {S,C,1Z}, is called weakly reversible if for any 
reaction v k — > v' k , there is a sequence of directed reactions beginning with v' k as a source complex 
and ending with v k as a product complex. That is, there exist complexes v\,...,v r such that 
v' k — > v\, v\ — \ v 2 , ■ ■ ■ , Vr — >■ v k & 1Z. A network is called reversible if v' k — > v k € 1Z whenever 

v k — > v' k e 1Z. 

Remark. The definition of a reversible network given in Definition 12.21 is distinct from the notion 
of a reversible stochastic process. However, in Section H31 we point out a connection between the 
two concepts for systems that are detailed balanced. 

To each reaction network, {S, C, 7Z}, there is a unique, directed graph constructed in the fol- 
lowing manner. The nodes of the graph are the complexes, C. A directed edge is then placed from 
complex v k to complex v' k if and only if v k — >■ v' k 6 1Z. Each connected component of the resulting 
graph is termed a linkage class of the graph. We denote the number of linkage classes by i. It is 
easy to see that a chemical reaction network is weakly reversible if and only if each of the linkage 
classes of its graph is strongly connected. 

Definition 2.3. S = span| Vfc ^^ e7 ^|{z/(, — u k } is the stoichiometric subspace of the network. For 
c G IR m we say c + S and (c + S) fl M> are the stoichiometric compatibility classes and positive 
stoichiometric compatibility classes of the network, respectively. Denote dim(S') = s. 

It is simple to show that for both stochastic and deterministic models, the state of the system 
remains within a single stoichiometric compatibility class for all time, assuming that one starts 
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in that class. This fact is important because it changes the types of questions that are reasonable 
to ask about a given system. For example, unless there is only one stoichiometric compatibility 
class, and so S = M m , the correct question is not whether there is a unique fixed point for a 
given deterministic system. Instead, the correct question is whether within each stoichiometric 
compatibility class there is a unique fixed point. Analogously, for stochastically modeled systems 
it is typically of interest to compute stationary distributions for each closed, irreducible subset of 
the state space (each contained within a stoichiometric compatibility class) with the precise subset 
being determined by initial conditions. 

The final definition of this section is that of the deficiency of a network (" lfT6lO . It is not a difficult 
exercise to show that the deficiency of a network is always greater than or equal to zero. 

Definition 2.4. The deficiency of a chemical reaction network, {S, C, 1Z}, is 5 = \C\ — £ — s, where 
\C\ is the number of complexes, £ is the number of linkage classes of the network graph, and s is 
the dimension of the stoichiometric subspace of the network. 

While the deficiency is, by definition, only a property of the network, we will see in Sections 
I3.2L |4l and [6] that a deficiency of zero has implications for the long-time dynamics of both deter- 
ministic and stochastic models of chemical reaction systems. 

3 Dynamical models 

The notion of a chemical reaction network is the same for both stochastic and deterministic sys- 
tems and the choice of whether to model the evolution of the state of the system stochastically or 
deterministically is made based upon the details of the specific chemical or biological problem at 
hand. Typically if the number of molecules is low, a stochastic model is used, and if the number of 
molecules is high, a deterministic model is used. For cases between the two extremes a diffusion 
approximation can be used or, for cases in which the system contains multiple scales, pieces of the 
reaction network can be modeled stochastically, while others can be modeled deterministically (or, 
more accurately, absolutely continuously with respect to time). See, for example, [8] and Section 

3.1 Stochastic models 

The simplest stochastic model for a chemical network {S, C, 1Z} treats the system as a continuous 
time Markov chain whose state X e Z> is a vector giving the number of molecules of each 
species present with each reaction modeled as a possible transition for the state. We assume a 
finite number of reactions. The model for the A;th reaction, v\. — > v' k , is determined by the vector 
of inputs, v k , specifying the number of molecules of each chemical species that are consumed in 
the reaction, the vector of outputs, v' k , specifying the number of molecules of each species that are 
created in the reaction, and a function of the state, Afc(X), that gives the rate at which the reaction 
occurs. Specifically, if the /cth reaction occurs at time t, the new state becomes 

X(t) = X(t-) + v' k - v k . 
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Let Rk(t) denote the number of times that the kth reaction occurs by time t. Then the state of the 
system at time t can be written as 



X(t) = X(0) + RkiWk ~ (3.1) 
k 

where we have summed over the reactions. The process R k is a counting process with intensity 
X k (X(t)) (called the propensity in the chemistry literature) and can be written as 

Rkit) = Y k Qf X k (X(s))ds^j , (3.2) 

where the Y k are independent, unit-rate Poisson processes ((321, iPPfll Ch. 1 1). Note that (13.11 ) and 
(13.21) give a system of stochastic equations that uniquely determines X up to sup{t : ^ k R k (t) < 
oo}. The generator for the Markov chain is the operator, A, defined by 

A f{x) = Afc(x)(/(x + i/ k - v k ) - f(x)), (3.3) 

k 

where / is any function defined on the state space. 

A commonly chosen form for the intensity functions X k is that of stochastic mass-action, which 
says that for x e Z> the rate of the kth reaction should be given by 

( rri \ / \ m I 

n J w = k * n j^yk^y, (3.4) 

for some constant K k , where we adopt the convention that 0! = 1. Note that the rate (|3.4I) is 
proportional to the number of distinct subsets of the molecules present that can form the inputs for 
the reaction. Intuitively, this assumption reflects the idea that the system is well-stirred in the sense 
that all molecules are equally likely to be at any location at any time. For concreteness, we will 
assume that the intensity functions satisfy (13.41) throughout most of the paper. In Section[6]we will 
generalize our results to systems with more general kinetics. 

A probability distribution {n(x)} is a stationary distribution for the chain if 

5>(x)A/(x) = 

x 

for a sufficiently large class of functions / or, taking f(y) = l x (y) and using equation (13.31) . if 

Y n ( x ~ v 'k + v k)X k {x -v' k + v k ) = ir(x) Y ^k(x) (3.5) 
k k 

for all x in the state space. If the network is weakly reversible, then the state space of the Markov 
chain is a union of closed, irreducible communicating classes. (This fact follows because if the 
Markov chain can proceed from state x to state y via a sequence of reactions, weak reversibil- 
ity of the network implies those reactions can be "undone" in reverse sequential order by another 
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sequence of reactions.) Also, each closed, irreducible communicating class is either finite or count- 
able. Therefore, if a stationary distribution with support on a single communicating class exists it 
is unique and 

lim P(X(t) = x | X(0) = y) = tt(z), 

t— >oo 

for all x, y in that communicating class. Thus, the stationary distribution gives the long-term 
behavior of the system. 

Solving equation (13.51) is in general a formidable task. However, in Section 0] we will do so if 
the network is weakly reversible, has a deficiency of zero, and if the rate functions Xk(x) satisfy 
mass-action kinetics, (13.41) . We will also show that the stationary distribution is of product form. 
More specifically, we will show that for each communicating class there exists a c G M> and a 
normalizing constant M > such that 

m m x 

7r(x) = Mj]vr J (x i ) = M]J%- 

• 1 ■ i 

1=1 t=l 

satisfies equation (|3.5I) . The q in the definition of Hi will be shown to be the zth component of an 
equilibrium value of the analogous deterministic system described in the next section. In Section 
[6] we will solve (13.51) for more general kinetics. 

3.2 Deterministic models and the deficiency zero theorem 

Under an appropriate scaling limit (see Section [47Tb the continuous time Markov chain (13.11 ), (13.21) , 
(13.41) becomes 

x (t) = x(0) + /^( s )) ds ) K - v k ) := x(0) + jf f(x(s))ds, (3.6) 

where the last equality is a definition and 

Jky-L) — rvk- L l x 2 x m 1 \ J -') 

where we use the convention 0° = 1. We say that the deterministic system (13.61) has mass-action 
kinetics if the rate functions ff. have the form (13 .7b ■ The proof of the following theorem by Feinberg 
can be found in ifToTl or [fT9ll . We note that the full statement of the deficiency zero theorem actually 
says more than what is given below and the interested reader is encouraged to see the original work. 

Theorem 3.1 (The Deficiency Zero Theorem). Consider a weakly reversible, deficiency zero chem- 
ical reaction network {S, C, 71} with dynamics given by (I3.6I) - (I3.7I) . Then for any choice of rate 
constants {nk}, within each positive stoichiometric compatibility class there is precisely one equi- 
librium value, and that equilibrium value is locally asymptotically stable relative to its compati- 
bility class. 

The dynamics of the system (I3.6I) - (I3.7I) take place in K> . However, to prove the deficiency 
zero theorem it turns out to be more appropriate to work in complex space, denoted M. c , which we 
will describe now. For any U C C let Uu : C — > {0, 1} denote the indicator function cuu(v k ) = 
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l{^ fe ec/}- Complex space is defined to be the vector space with basis {u Uk \ v k G C}, where we have 
denoted uJ{ Uk } by uo Uk . 

If u is a vector with nonnegative integer components and w is a vector with nonnegative real 
components, then let u\ = Yli u ^- an ^ w u = Yli w i^ where we interpret 0° = 1 and 0! = 1. Let 

* : R rn ->■ R c and A R : R c -»■ M c be defined by: 

where the subscript k of A K denotes the choice of rate constants for the system. Let Y : R c — > R m 
be the linear map whose action on the basis elements {uJ Uk } is defined by Y(u Vk ) = v k . Then 
equations (13 .61) - (13 .7 1) can be written as the coupled set of ordinary differential equations 

x(t) = f(x(t))=Y(A K (V(x(t)))). 

Therefore, in order to show that a value c is an equilibrium of the system, it is sufficient to show 
that A K ($l(c)) = 0, which is an explicit system of equations for c. In particular, A k (^/(c)) = if 
and only if for each z G C 

{k:u' k =z} {k:u k =z} 

where the sum on the left is over reactions for which z is the product complex and the sum on the 
right is over reactions for which z is the source complex. 

The following has been shown in ll28l and lfT6l (see also |]25l0 . 

Theorem 3.2. Let {S, C, 71} be a chemical reaction network with dynamics given by (I3.6I) - (I3.7I) 

for some choice of rate constants, {^k}- Suppose there exists a c G R™ for which A K (^(c)) = 0, 
then the following hold: 

1. The network is weakly reversible. 

2. Every equilibrium point with strictly positive components, x G M> with fix) = 0, satisfies 
A K {V{x)) =0. 

3. If Z = {x G M™ | f{x) = 0}, then \nZ := {y G R m | 3 x G Z and y { = m^)} is 
a coset of S*- 1 , the perpendicular complement of S. That is, there is a k G R m such that 
In Z = {w G R m | w = k + ufor some u G 5'" L }. 

4. There is one, and only one, equilibrium point in each positive stoichiometric compatibility 
class. 

5. Each equilibrium point of a positive stoichiometric compatibility class is locally asymptoti- 
cally stable relative to its stoichiometric compatibility class. 
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Thus, after a choice of rate constants has been made, the conclusions of the deficiency zero 
theorem pertaining to the existence and asymptotic stability of equilibria (points 4. and 5. of The- 
orem |3.2|) hold so long as there exists at least one c e M.™ such that A K (^f(c)) = 0. The condition 
that the system has a deficiency of zero only plays a role in showing that there does exist such a 
c G K™ . A proof of the following can be found in |fT6ll, ifTTll, or lfl9l. 

Theorem 3.3. Let {S, C, 1Z} be a chemical reaction network with dynamics given by (|3.6I) - (|3.7I) 

for some choice of rate constants, {ku}- If the network has a deficiency of zero, then there exists a 
c G M> such that A K (^(c)) = if and only if the network is weakly reversible. 

A chemical reaction network with deterministic mass-action kinetics (and a choice of rate con- 
stants) that admits a c for which A K (^(c)) = is called complex balanced in the literature. The 
second conclusion of Theorem 13 .21 demonstrates why this notation is appropriate. The equivalent 
representation given by (13.81) shows the origin of this terminology. The surprising aspect of the 
deficiency zero theorem is that it gives simple and checkable sufficient conditions on the network 
structure alone that guarantee that a system is complex balanced for any choice of rate constants. 
We will see in the following sections that the main results of this paper have the same property: 
product-form stationary distributions exist for all stochastic systems that are complex balanced 
when viewed as deterministic systems, and 5 = is a sufficient condition to guarantee this for 
weakly reversible networks. 



4 Main result for mass-action systems 

The collection of stationary distributions for a countable state space Markov chain is convex. The 
extremal distributions correspond to the closed, irreducible subsets of the state space; that is, every 
stationary distribution can be written as 

7r = ^a r vr r , (4.1) 
r 

where a r > 0, J2r a r = 1> and the sums are over the closed, irreducible subsets T of the state 
space. Here 7r r is the unique stationary distribution satisfying 7r r (r) = 1. 

We now state and prove our main result for systems with mass-action kinetics. 

Theorem 4.1. Let {S,C,TZ} be a chemical reaction network and let {nk} be a choice of rate 
constants. Suppose that, modeled deterministically, the system is complex balanced with complex 
balanced equilibrium c G M> . Then the stochastically modeled system with intensities (13.41) has a 
stationary distribution consisting of the product ofPoisson distributions, 

m Xi 

i=l Xi - 

If Z> is irreducible, then (|4.2I) is the unique stationary distribution, whereas if Z> is not irre- 
ducible then the of equation (14.11) are given by the product-form stationary distributions 

m x . 

7r r (x) = M r xer, 

1=1 

and 7Tr(x) = otherwise, where Mr is a positive normalizing constant. 
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Proof. Let tt satisfy (14.21 ) where c G M> satisfies A K (^(c)) = 0. We will show that n is stationary 
by verifying that equation (13.51 ) holds for all x G Z> . Plugging 7r and (13.41 ) into equation (13.51) and 
simplifying yields 

e /t ^" < (x _^) ! n = e Kfc (7T^)T n i^^}- ( 4 -3) 

Equation (14.31 ) will be satisfied if for each complex zgC, 

j m ^ m 

E ^^T^^iII 1 ^^ E ^T^yll 1 ^}' (4-4) 

{fc:^= 2 } V ' '' 1=1 {k:u k =z} V ' ' 1=1 

where the sum on the left is over reactions for which z is the product complex and the sum on 
the right is over reactions for which z is the source complex. The complex z is fixed in the above 
equation, and so (14.41) is equivalent to (13.81) . which is equivalent to A K (^(c)) = 0. 

To complete the proof, one need only observe that the normalized restriction of n to any closed, 
irreducible subset T must also be a stationary distribution. □ 

The following theorem gives simple and checkable conditions that guarantee the existence of a 
product-form stationary distribution of the form (|4.2I) . 

Theorem 4.2. Let {S, C, 71} be a chemical reaction network that has a deficiency of zero and is 
weakly reversible. Then for any choice of rate constants {n k } the stochastically modeled system 
with intensities (13.41) has a stationary distribution consisting of the product of Poisson distributions, 



1=1 



where c is an equilibrium value for the deterministic system (|3.6I) - (|3.7I) . which is guaranteed to 
exist and be complex balanced by Theorems \3.1^3.3\ If Z> is irreducible, then tt is the unique 



stationary distribution, whereas j/Z> is not irreducible then the of equation (|4.1I) are given by 
the product-form stationary distributions 



m Xi 



7rr(x) = M r xeT, 



-. x% ■ 
1=1 



and 7r r (x) = otherwise, where M r is a positive normalizing constant. 

Proof. This is a direct result of Theorems 13.31 and 14.11 □ 

We remark that Theorems 14.11 and 14.21 give sufficient conditions under which Z> being irre- 
ducible guarantees that when in distributional equilibrium the species numbers: (a) are independent 
and (b) have Poisson distributions. We return to this point in Examples 15 . 21 and | 5 . 3 1 
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4.1 The classical scaling 

Defining \j/k\ = X^i u %k an d letting V be a scaling parameter usually taken to be the volume of 
the system times Avogadro's number, it is reasonable to scale the rate constants of the stochastic 
model with the volume like 

«* = yfip;, (4-5) 

for some k& > 0. This follows by considering the probability of a particular set of \i/k\ molecules 
finding each other in a volume proportional to V in a time interval [t, t + At). In this case, the 
intensity functions become 

Since V is the volume times Avogadro's number and x gives the number of molecules of each 
species present, c = V~ l x gives the concentrations in moles per unit volume. With this scaling 
and a large volume limit 

\l{x) « V7c* J] ^ = Vkk ° Uk = V ^(c)- (4.7) 

i 

Since the law of large numbers for the Poisson process implies V~ l Y k (Nu) ~ u, (13.21 ) and (14.71 ), 
together with the assumption that X(0) = VC(0) for some C(0) G M> , imply 

C(t) = V- l X(t) « C(0) + V / « fc C(s)"*tfe K - z/ fe ), 

which in the large volume limit gives the classical deterministic law of mass action detailed in Sec- 
tion [372]. For a precise formulation of the above scaling argument, termed the "classical scaling," 

Because the above scaling is the natural relationship between the stochastic and deterministic 
models of chemical reaction networks, we expect to be able to generalize Theorem 14.11 to this 
setting. 

Theorem 4.3. Let {S, C, 71} be a chemical reaction network. Suppose that, modeled deterministi- 
cally with rate constants the system is complex balanced with complex balanced equilibrium 
c G R> . For some V > 0, let {n^} be related to {k/-} via (14.51) . Then the stochastically modeled 
system with intensities (13.41 ) and rate constants {n k } has a stationary distribution consisting of the 
product of Poisson distributions, 



T ■ I 

i=i l " 



If Z> is irreducible, then (|4.2I) is the unique stationary distribution, whereas if Z> is not irre- 
ducible then the 7Tr of equation (14.11) are given by the product-form stationary distributions 

vr r (x) = M r n^T- ? ^T, 
and 7Tr(x) = otherwise, where Mr is a positive normalizing constant. 
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Proof. The proof is similar to before, and now consists of making sure the V's cancel in an appro- 
priate manner. The details are omitted. □ 



We see that Theorem 14 . 1 1 folio w s from Theorem 14.3 1 by taking V — 1. Theorem l4.2l generalizes 
in the obvious way. 

4.2 Reversibility and detail balance 

An equilibrium value, c E K.> , for a reversible, in the sense of Definition I2.2L chemical reaction 
network with deterministic mass-action kinetics is called detailed balanced if for each pair of 
reversible reactions, v k +± u' k , we have 

n k c Vk = 4c^, (4.8) 

where K k , n' k are the rate constants for the reactions v k — > u' k , v' k — V v k , respectively. In [fT8l , page 
1820, Feinberg shows that if one positive equilibrium is detailed balanced then they all are; a re- 
sult similar to the second conclusion of Theorem 13 .21 for complex balanced systems. A reversible 
chemical reaction system with deterministic mass action kinetics is therefore called detailed bal- 
anced if it admits one detailed balanced equilibrium. It is immediate that any system that is detailed 
balanced is also complex balanced. The fact that a product-form stationary distribution of the form 
(14.21) exists for the stochastic systems whose deterministic analogs are detailed balanced is well- 
known. See, for example, 11381 . Theorems 14. II and 14.21 can therefore be viewed as an extension 
of that result. However, more can be said in the case when the deterministic system is detailed 
balanced, and which we include here for completeness (no originality is being claimed). 

As mentioned in the remark following Definition 12.21 the term "reversible" has a meaning in 
the context of stochastic processes that differs from that of Definition 12 .21 Before defining this, we 
need the concept of a transition rate. For any continuous time Markov chain with state space T, the 
transition rate from x E T to y E T (with x ^ y) is a non-negative number a(x, y) satisfying 

P(X(t + At) = y | X{t) = x)= a(x, y)At + o(At). 

Thus, in the context of this paper, if y = x + v' k — v k for some k, then a(x, y) = X k (x), and zero 
otherwise. 

Definition 4.4. A continuous time Markov chain X(t) with transition rates a(x, y) is reversible 
with respect to the distribution n if for all x, y in the state space T 

n(x)a(x, y) = n(y)a(y, x). (4.9) 

It is simple to see (by summing both sides of (14.91) with respect to y over T), that n must be 
a stationary distribution for the process. A stationary distribution satisfying (14.91) is even called 
detailed balanced in the probability literature. The following is proved in [|38l . Chapter 7. 

Theorem 4.5. Let {S, C, 71} be a reversible 2 chemical reaction network with rate constants {n k }. 
Then the deterministically modeled system with mass-action kinetics has a detailed balanced equi- 
librium if and only if the stochastically modeled system with intensities (13.41) is reversible with 
respect to its stationary distribution} 
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Succinctly, this theorem says that reversibility and detailed balanced in the deterministic setting 
is equivalent to reversible (and, hence, detailed balanced) in the stochastic setting. 



4.3 Non-uniqueness of c 

For stochastically modeled chemical reaction systems any irreducible subset of the state space, 
T, is contained within (y + S) D Z> for some y E M> . Therefore, each T is associated with 
a stoichiometric compatibility class. For weakly reversible systems with a deficiency of zero, 
Theorems 13 .21 and [3731 guarantee that each such stoichiometric compatibility class has an associated 
equilibrium value for which A K (^(c)) = 0. However, neither Theorem 14.11 nor Theorem 14.21 
makes the requirement that the equilibrium value used in the product-form stationary measure 
7Tr(-) be contained within the stoichiometric compatibility class associated with T. Therefore 
we see that one such c can be used to construct a product-form stationary distribution for every 
closed, irreducible subset. Conversely, for a given irreducible subset T any positive equilibrium 
value of the system (I3.6I) - (I3.7I) can be used to construct 7r r (-). This fact seems to be contrary 
to the uniqueness of the stationary distribution; however, it can be understood through the third 
conclusion of Theorem 13 .21 as follows. 

Let T be a closed, irreducible subset of the state space with associated positive stoichiometric 
compatibility class (y + S) n Z^ , and let c h c 2 E M> be such that A K (^(ci)) = A K (^(c 2 )) = 0. 
For i E {1,2} and x E T, let %i(x) = MiC* / x\, where Mi and M 2 are normalizing constants. Then 
for each x E T 

tti(x) _ Micf x\ M x c\ 

7T 2 (x) ~ X\ M 2 C% ~ M2 Cf ' 

For any vector u, we define (In (it)), = ln(wj). Then for x E Y C y + S 

£l _ e 2>(lnci-lne 2 ) _ ^-(lnci-lnca) _ £l (A -^Q) 

where the second equality follows from the third conclusion of Theorem l3.2l Therefore, 

7Ti(x) _ Mx C\ 

Mxj ~ W 2 4' 

Finally, 



(4.11) 



1= (m 1 J2cVx\)/(m 2 J2c x 2 /xi\ 
\ xer J \ x£F J 



7T 2 (x) 

where the second equality follows from equation (|4.10l) and the third equality follows from equa- 
tion (14.1 II) . We therefore see that the stationary measure is independent of the choice of c, as 
expected. 



2 In the sense of Definition ^. 21 
3 in the sense of Definition E4] 
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5 Examples 



Our first example points out that the existence of a product-form stationary distribution for the 
closed, irreducible subsets of the state space does not necessarily imply independence of the species 
numbers. 

Example 5.1. (Non-independence of species numbers) Consider the simple reversible system 

kx 

Si ^ S 2 , 
k 2 

where k\ and k 2 are nonzero rate constants. We suppose that Xi(Q) + X 2 (0) = N, and so X%(t) + 
X 2 (t) = N for all t. This system has two complexes, one linkage class, and the dimension of 
the stoichiometric compatibility class is one. Therefore it has a deficiency of zero. Since it is also 
weakly reversible, our results hold. An equilibrium to the system that satisfies the complex balance 
equation is 

k 2 ki 



ki + k 2 k x + k 2 
and the product-form stationary distribution for the system is 



^1 ^2 
n(x) = M- 



.C\ c 2 



X\\ x 2 \ 

where M > is a normalizing constant. Using that X\(t) + X 2 (t) = N for all t yields 

Xi N-Xi yr 

K ' x x \[X-x x )\ x x \[X-x x )\ 1 v 1 
for < X\ < N. After setting M = N\, we see that X\ is binomially distributed. Similarly, 

^ 2/ 



for < x 2 < N. Therefore, we trivially have that P(X t = N) = cf and P(X 2 = N) = eg, but 
P(Xx = N, X 2 = N) = ^ cf eg, and so X x and X 2 are not independent. 

Remark. The conclusion of the previous example, that independence does not follow from the exis- 
tence of a product-form stationary distribution, extends trivially to any network with a conservation 
relation among the species. 

Example 5.2. (First order reaction networks) The results presented below for first order reaction 
networks are known in both the queueing theory and mathematical chemistry literature. See, for 
example, Il3~0ll and l|20l . We present them here to point out how they follow directly from Theorem 

We begin by defining |i>| = X^ 1 "* f° r an Y vector v E IR> . We say a reaction network is a 
first order reaction network if \v^.\ G {0, 1} for each complex v\. E C. Therefore, a network is first 
order if each entry of the v k are zeros or ones, and at most one entry can be a one. It is simple to 
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show that first order reaction networks necessarily have a deficiency of zero. Therefore, the results 
of this paper are applicable to all first order reaction networks that are weakly reversible. Consider 
such a reaction network with only one linkage class (for if there is more than one linkage class we 
may consider the different linkage classes as distinct networks). We say that the network is open if 
there is at least one reaction, — > u' k , for which = 0. Hence, by weak reversibility, there must 
also be a reaction for which v' k = 0. If no such reaction exists, we say the network is closed. If the 



network is open we see that S 
stationary distribution is 



T = Z> is irreducible, and so by Theorem 14.21 the unique 



7T (X) 



n 



x G Z 



>0; 



where c G M> is the complexed balanced equilibrium of the associated (linear) deterministic 
system. Therefore, when in distributional equilibrium, the species numbers are independent and 
have Poisson distributions. Note that neither the independence nor the Poisson distribution resulted 
from the fact that the system under consideration was a first order system. Instead both facts 
followed from T being all of Z> . 

In the case of a closed, weakly reversible, single linkage class, first order reaction network, it 
is easy to see that there is a unique conservation relation Xi(t) + ■ — h X m {t) = N, for some N. 
Thus, in distributional equilibrium X{t) has a multinomial distribution. That is for any x G Z' 



>o 



satisfying x 1 + x 2 + ■ ■ ■ + x Ti 



ir(x) 



N 



N 



Xt,X 2 , 



x r . 



Xi\ ■■■Xr. 



(5.1) 



where c G M™ is the equilibrium of the associated deterministic system normalized so that £\ q = 
1. As in the case of the open network, we note that the form of the equilibrium distribution does 
not follow from the fact that the network only has first order reactions. Instead (15.11) follows from 
the structure of the closed, irreducible communicating classes. 

Example 5.3. (Enzyme kinetics I) Consider the possible model of enzyme kinetics given by 

E + S^ES^E + P, E ^ +±S, (5.2) 

where E represents an enzyme, S represents a substrate, ES represents an enzyme- substrate com- 
plex, P represents a product, and some choice of rate constants has been made. We note that both 
E and S are being allowed to enter and leave the system. 



The network (15.21) is reversible and has six complexes and two linkage classes. The dimension 
of the stoichiometric subspace is readily checked to be four, and so the network has a deficiency of 
zero. Theorem 14 . 21 applies and so the stochastically modeled system has a product-form stationary 
distribution of the form (14.21) . Ordering the species as X\ = E, X 2 = S, X 3 = ES, and X 4 = P, 
the reaction vectors for this system include 



1 








1 






-1 




1 


-1 







1 


1 


-1 







1 
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We therefore see that T = Z> is the unique closed, irreducible communicating class of the 
stochastically modeled system and Theorem 14.21 tells us that in distributional equilibrium the 
species numbers are independent and have Poisson distributions with parameters q, which are 
the complex balanced equilibrium values of the analogous deterministically modeled system. 

Example 5.4. (Enzyme kinetics II) Consider the possible model for enzyme kinetics given by 

ki k2 ks 

E + S^ES^E + P, ^ £, (5.3) 

k-l fc-2 

where the species E, S, ES, and P are as in Example l5.3l We are now allowing only the enzyme 
E to enter and leave the system. The network is reversible, there are five complexes, two linkage 
classes, and the dimension of the stoichiometric compatibility class is three. Therefore, Theorem 
|4.2| implies that the stochastically modeled system has a product-form stationary distribution of the 
form (14.21) . The only conserved quantity of the system is S + ES + P, and so X 2 [t) + X 3 (t) + 
X 4 (t) = iV for some N > and all t. Therefore, after solving for the normalizing constant, we 
have that for any x G Z> satisfying x 2 + x 3 + x 4 = N 

r Xl /VI r Xl f N \ 

_/„\ — -CI 1 X2 Xs Xa — Ci If 1 Xo Xs Xa 



x\\ x^-Xy-x^. x\\ \x 2 , X3, X4 

where c = (k 3 /k_ 3 , c 2 , C3, C4) has been chosen so that 02 + 03 + 04 = 1. Thus, when the stochasti- 
cally modeled system is in distributional equilibrium we have that: (a) E has a Poisson distribution 
with parameter fc 3 /fc_ 3 , (b) S, ES, and P are multinomially distributed, and (c) E is independent 
from S, ES, and P. 



5.1 The multiscale nature of reaction networks 

Within a cell, some chemical species may be present in much greater abundance than others. In 
addition, the rate constants Kk may vary over several orders of magnitude. Consequently, the 
scaling limit that gives the classical deterministic law of mass action detailed in Section |4~T1 may 
not be appropriate, and a different approach to deriving a scaling limit approximation for the basic 
Markov chain model must be considered. As a consequence of the multiple scales in a network 
model, it may be possible to separate the network into subnetworks of species and reactions, each 
dominated by a time scale of a specific magnitude. Within each subnetwork, the graph structure and 
stoichiometric properties may determine properties of the asymptotic solutions of the subnetwork. 

Example 5.5. Consider the reaction network 

S + E^C^P + Eu E 1 ^A + E 2 , 0^£ 2 , 

where — ► E 2 and E 2 — > represent production and degradation of E 2 , respectively, S is a 
substrate being converted to a product P, E\ and E 2 are enzymes, and A is a substrate that reacts 
with E 2 allosterically to transform it into an active form. 

We suppose that {%) the enzymes E\, E 2 and the substrate A are in relatively low abundances, 
(ii) the substrate S has a large abundance of 0(V), and (Hi) the reaction rates are also of the order 
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0{V). We change notation slightly and denote the number of molecules of species A at time t as 
X A (t), and similarly for the other species. Further, we denote Xg(t) /V = Zg(t). Combined with 
the conservation relation X E + X^, + X A = M G Z >0 , the scaled equations for the stochastic 
model are 

Z v s (t) = Z v s (V)-V- 1 Y l {V [ K l Z v s {s)Xl i {s)ds) + V- l Y 2 {V [ K 2 X%(s)ds) 

Jo Jo 

Xl(t)=Xl(0)-Y 1 (V f K l Z v s {s)X v Ei {s)ds) + Y 2 {V [ n 2 X^{s)ds) 

Jo Jo 

+ Y 3 (V [ K 3 X^(s)ds)-Y 4 (V [ K4 X^(s)ds) +Y 5 (V f K b X v A {s)X v E2 {s)ds) 
Jo Jo Jo 

X v A {t)=X v A {Q) + Y A {V [ K A X v Ei (s)ds)-Y,(V [ ^X v A (s)X v E2 (s)ds) 

Jo Jo 

X^ 2 (t)=Xl 2 (0) + Y 6 (V^t) + Y 4 (V [ K A Xl 1 {s)ds)-Y 5 {V I ^X v A {s)X v E2 {s)ds) 

Jo Jo 

-Y 7 (V [ K 7 X^{s)ds), 
Jo 

where the Yi are unit-rate Poisson processes. The first equation satisfies 

pt poo pt poo 

Z%(t) = Zj(0) - V^YiiV / KtZ^s) / x^(dx)ds) + V~ l Y 2 (V k 2 xr%(dx)ds), 

Jo J -oo Jo J-oo 

where fi^(A) = I^ x v ( s ) e A} an d vY(A) = I{x v (s)eA} are the respective occupation measures. 
Using methods from stochastic averaging (see, for example, CUES), as V — > oo the fast system is 
"averaged out:" 

pt poo pt poo 

Zs(t) = Zs(0) — / KiZs(s) / x[i s (dx)ds + / k 2 I xr] s {dx)ds 1 (5.4) 

Jo J-oo Jo J-oo 

where fi s and r] s are the stationary distributions of X El and X c , respectively, of the fast subsystem 
with Zs(s) held constant (assuming a stationary distribution exists). This reduced network (i.e. the 
fast subsystem) is 

A + E 2 ^±E 1 i± C, 0^£ 2 . (5.5) 

Setting z = Z s (s) we have the following equilibrium relations for the moments of the above 
network 

k 4 E[X Ei ] - k 5 E[X a X E2 ] = 
-{k±z + k a )E[X Ei ] + («2 + k 3 )M[X c ] + k 5 E[X a X E2 ] = 
k 6 + n 4 E[X El ] - k 5 E[X a X E2 ] - k 7 E[X E2 ] = 
E[X El ] + E[X C ] + E[X A ] = M. 

E[XeJ and E[Xc], which are both functions of z and needed in equation (15.41) . can not be explic- 
itly solved for via the above equations without extra tools as (15.61) is a system of four equations 
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with five unknowns. This situation arises frequently as it stems from the nonlinearity of the sys- 
tem. However, the network (15.51 ) consists of five complexes, two connected components, and the 
dimension of its stoichiometric subspace is three. Therefore, its deficiency is zero. As it is clearly 
weakly reversible, Theorem 14.11 applies and, due to the product form of the distribution and the 
unboundedness of the support of Xe 2 , it is easy to argue that Xe 2 is independent of Xa, Xe x , and 
X c when in equilibrium. Thus, EfX^X^] = E[X^]E[X S2 ] and the first moments can be solved 
for as functions of Z s (s). After solving and inserting these moments, (15.41 ) becomes 

k 1 k 3 k 5 k 6 MZ s (s) ^ 

which is Michaelis-Menten kinetics. 



Z s (t) = Z s (0) - 



6 More general kinetics 

In this section we extend our results to systems with more general kinetics than stochastic mass 
action. The generalizations we make are more or less standard for the types of results presented 
in this paper (see, for example, lf30ll . Section 8.5, OH , Chapter 9). What is surprising, however, is 
that the conditions of the deficiency zero theorem of Feinberg (which are conditions on mass-action 
deterministic systems) are also sufficient to guarantee the existence of stationary distributions of 
stochastically modeled systems even when the intensity functions are not given by (I3.4I ). It is 
interesting to note that the generalizations made here for the stochastic deficiency zero Theorem 
I4.2l are similar to those made in ||37ll . which generalized Feinberg 's deficiency zero Theorem l3.1l in 
the deterministic setting. 

Suppose that the intensity functions of a stochastically modeled system are given by 

m v ik -l m 

Afc(x) K/, J) ] [ 9i(xi - j) = n k Y[ 9i{xi)ei(xi - l)9i(xi - [y ik - 1)), (6.1) 

i=l j=0 i=l 

where the K k are positive constants, 0j : Z — > M>o, Oi(x) = if x < 0, and we use the convention 
that ri7=o a i = 1 f° r an y { a A- Note that the final condition allows us to drop the indicator 
functions of (|3.4I) . As pointed out in 11301 . the function 8{ should be thought of as the "rate of 
association" of the ith species. We give a few interesting choices for 6>j. If 6>j(xj) = Xj for Xi > 0, 
then (16.11) is stochastic mass-action kinetics. However, if for Xi > 

B^) = t~~~~ > (6-2) 

f^i \ Xi 

for some positive constants ki and Vi, then the system has a type of stochastic Michaelis-Menten 
kinetics ([29], Chapter 1). Finally, if \u k \ 6 {0, 1} and 6i(xi) = min{rij, Xi] for Xi > 0, then the 
dynamical system models an M/M/n queueing network in which the ith species (and in this case 
complex) represents the queue length of the ith queue, which has n; servers who work on a first 
come, first serve basis. 

The main restriction imposed by (16.11 ) is that for any reaction for which the ith species appears 
in the source complex, the rate of that reaction must depend upon X { via 9i(Xi) only. Therefore, 
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if, say, the ith species is governed by the kinetics (16.21 ), then the constants ki and Vi must be the 
same for each intensity which depends upon X { (although the v,i may be incorporated into the rate 
constants K k , and so the real restriction is on the constant ki). However, systems with intensities 
given by (16.11) are quite general in that different kinetics can be incorporated into the same model 
through the functions 6>j. For example, if in a certain system species Si is modeled to be governed 
by Michaelis-Menten kinetics (16.21) and species S 2 is modeled to be governed by mass-action 
kinetics, then the reaction Si + 5 2 — > v' k would have intensity 

ViXi 

X k {X) = K k - ; X 2 , 

ki + xi 

for some constant n k . 

In following we use the convention that Ylj=i a j = 1 f° r an Y cn °i ce of {<%}• 

Theorem 6.1. Let {S, C, 71} be a stochastically modeled chemical reaction network with intensity 
functions (16.11) . Suppose that the associated mass-action deterministic system with rate constants 
{n k } has a complex balanced equilibrium c G R> . Then the stochastically modeled system admits 
the stationary distribution 

m Xi 

k(x) = M\[ 1 xGZ™, (6.3) 

i= i llj=i t wJ 

where M > is a normalizing constant, provided that (|6.3I) is summable. If Z> is irreducible, 
then (|6.3I) is the unique stationary distribution, whereas if Z> is not irreducible then the 7ir of 
equation (14.11) are given by the product-form stationary distributions 

m x . 

7r r (x) = M r n U x^ ei{j y ^T, (6.4) 



and 7ir(x) = otherwise, where Mr > is a normalizing constant, provided that (|6.4I) is 
summable. 

Proof. The proof consists of plugging (16.31) and (16.11 ) into equation (13.51 ) and verifying that c being 
a complex balanced equilibrium is sufficient. The details are similar to before and so are omitted. 

□ 

Remark. We simply remark that just as Theorem 14.21 followed directly from Theorem 14. 1[ the re- 
sults of Theorem 16. 1 1 hold, independent of the choice of rate constants n k , so long as the associated 
network is weakly reversible and has a deficiency of zero. 

Example 6.2. Consider a network, {S,C, 1Z}, that is weakly reversible and has a deficiency of 
zero. Suppose we have modeled the dynamics stochastically with intensity functions given by 
(16.11) with each 9i given via (|6.2I) for some choice of v j > and ki a nonnegative integer. That is, 
we consider a system endowed with stochastic Michaelis-Menten kinetics. Then, 



n*w=n^=<7 (T) 

3=1 3=1 J V 7 
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Thus, our candidate for a stationary distribution is 



m x m /r \ / \ 

^nn^-n(^)© 

Noting that 




= 0(x l ° i ), Xi -> oo, 



we see that tt(x) given by (16.51) is summable if Cj < v; L for each species Si whose possible abun- 
dances are unbounded. In this case, (|6.5I) is indeed a stationary distribution for the system. We 
note that the condition q < Vi for each species Si is both necessary and sufficient to guarantee 
summability if Z> is irreducible, as in such a situation the species numbers are independent. 

Example 6.3. In [1331 , Levine and Hwa computed and analyzed the stationary distributions of 
different stochastically modeled chemical reaction systems with Michaelis-Menten kinetics (16.21) . 
The models they considered included among others: directed pathways (0 — > Si — > S2 — > ■ ■ ■ — > 
Sl 0)» reversible pathways (0 — > Si S 2 • • • Sl — > 0), pathways with dilution of 
intermediates (Si — > 0), and cyclic pathways (S^ — >■ Si). Each of the models considered in [[331 
is biologically motivated and has a first order reaction network (\is k \ G {0, 1}, see Example 15 .21 ), 
which guarantees that they have a deficiency of zero. Further, the networks of the models consid- 
ered are weakly reversible; therefore, the results of the current paper, and in particular Theorem 
16.11 and the remark that follows, apply so long as the restrictions discussed in the paragraph pre- 
ceding Theorem 16. II are met. While these restriction are not always met (for example, dilution is 
typically modeled with a linear intensity function and there is no reason for the ki of a forward 
and a backward reaction for a species 5^ in a reversible pathway to be the same), they found that 
the stationary distributions for these models are either of product form (when the restrictions are 
met) or near product form (when the restrictions are not met). Further, because Z> is irreducible 
in each of these models, the product form of the distribution implies that the species numbers are 
independent. It is then postulated that the independence of the species numbers could play an 
important, beneficial, biological role (see ll35l for details). Similar to the conclusions we drew in 
Example 15 .21 Theorem 16. II and the remark that follows point out how the models analyzed in ll35l 
are actually special cases of a quite general family of systems that have both the product form and 
independence properties, and that these properties may be more widespread, and taken advantage 
of by living organisms, than previously thought. 

We return to the result of Example [6]2] pertaining to the summability of (16.51) and show that this 
can be generalized in the following manner. 

Theorem 6.4. Suppose that for some closed, irreducible Y C Z> , 7r r : Y — > IR> satisfies 

for some c G M™ and M > 0, where Qi : Z> — > M>o/or each i. Then 7Tp(x) is summable if for 
each i for which sup{xj | x G Y} = 00 we have thatOi(j) > Ci + e for some e > and j sufficiently 
large. 
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Proof. The conditions of the theorem immediately imply that there are positive constants C and p 
for which 7r r (x) < Ce~ p ^, for all x E T, which implies that %r(x) is summable. □ 

It is tempting to believe that the conditions of Theorem 16 .41 are in fact necessary, as in the case 
when Z> is irreducible. The following simple example shows this not to be the case. 

Example 6.5. Consider the reaction system with network 

<=t Sx + S 2 , 

where the rate of the reaction — > Si + S 2 is Xi(x) = 1, and the rate of the reaction Si + S 2 — > 

is X 2 (x) = 1 x 9i(xi)9 2 (x 2 ), where 

Assume further that Xl(0) = X 2 (0). For the more physically minded readers, we note that this 
model could describe a reaction system for which there is a chemical complex C = SiS 2 that 
sporadically breaks into its chemical constituents, which may then re-form. The complex C may 
be present in such high numbers relative to free S\ and £2 that we choose to model it as fixed, 
which leads to the above reaction network. 

We note that in this case, the reaction rates {tz k } for the corresponding mass-action deter- 
ministic system are both equal to one, and so an equilibrium value guaranteed to exist for the 
deterministically modeled system by the deficiency zero theorem is c = (1, 1). This system 
does not satisfy the assumptions of Theorem 16.41 because both Xi and X 2 are unbounded and 
linx^oo 9 2 (j) = 1/2 < 1 = c 2 . However, for any x eT = {x G Z> : xi = x 2 }, 

l + xi\ fl\ xi fl + x 2 \ ( 1 V' 2 /l + xi\ 2 f^ Xl 



x i J \& J V x 2 / vli/ 2 )/ V x i 

which is summable over Y. 

For the most general kinetics handled in this paper, we let the intensity functions of a stochas- 
tically modeled system be given by 

rn 



Xk ^ = Kk e( x -u k ) n k^e.}, (6-6) 



6{x 

where the K k are positive constants, and 9 : Z m — > R >0 . Note that if 

m Xi 

W = 1111^0'). 

i=l j=l 

for some functions 6>, ; , then (16.61) is equivalent to (16.11) . and so the following theorem implies The- 
orem |6.1[ It's proof is similar to the previous theorems and so is omitted. 
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Theorem 6.6. Let {S, C, 71} be a stochastically modeled chemical reaction network with intensity 
functions (I6.6I ). Suppose that the associated mass-action deterministic system with rate constants 
{ftfc} has a complex balanced equilibrium c G K> - Then the stochastically modeled system admits 
the stationary distribution 

ir( x ) = M-—\[c*\ (6.7) 

where M > is a normalizing constant, provided that (16.71) is summable. If Z> is irreducible, 
then (|6.7I) is the unique stationary distribution, whereas if Z™ is not irreducible then the 7r r of 
equation (14.11 ) are given by the product-form stationary distributions 

7Tr (x) = M T — -He?, xeT, (6.8) 

and vr r (x) = otherwise, where M r > is a normalizing constant, provided that (|6.8I) is 
summable. 

Remark. Similar to the remark following Theorem l6.ll we point out that the results of Theorem l6.6l 
hold, independent of the choice of rate constants K k , so long as the associated network is weakly 
reversible and has a deficiency of zero. 
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